#include "nfem_sf.h"

double jacoby2(double x00, double x01, double x10, double x11)
{
	return x00*x11-x01*x10;
}

double jacoby3(double x00, double x01, double x02, 
	double x10, double x11, double x12, 
	double x20, double x21, double x22)
{
	return x00*x11*x22+x01*x12*x20+x10*x21*x02-x02*x11*x20-x01*x10*x22-x00*x12*x21;
}

bool nfem_sf::linear::line_element_is_valid(double x, double x1, double x2)
{
	if (x1 >= x2)
	{
		throw std::invalid_argument("nfem_sf::linear::line_element_is_valid(...) got invalid definition range.");
		return false;
	}

	if (x < x1 || x > x2)
	{
		throw std::invalid_argument("nfem_sf::linear::line_element_is_valid(...) got an invalid coordinate.");
		return false;
	}

	return true;
}

double nfem_sf::linear::line(double x, double x1, double x2, unsigned int idx, valtype_e vt)
{
	switch (10*idx + vt)
	{
		case 0:
			return (x2 - x)/(x2 - x1);
		case 1:
			return -1.0/(x2 - x1);
		case 10:
			return (x - x1)/(x2 - x1);
		case 11:
			return 1.0/(x2 - x1);
		default:
			throw std::invalid_argument("nfem_sf::linear::line(...) got an invalid calculating type.");
			return 0.0;
	}
}

bool nfem_sf::linear::quad_element_is_valid(double x, double y, double x1, double x2, double y1, double y2)
{
	if (x1 >= x2 || y1 >= y2)
	{
		throw std::invalid_argument("nfem_sf::linear::quad_element_is_valid(...) got invalid definition range.");
		return false;
	}

	if (x < x1 || x > x2 || y < y1 || y > y2)
	{
		throw std::invalid_argument("nfem_sf::linear::quad_element_is_valid(...) got an invalid coordinate.");
		return false;
	}

	return true;
}

double nfem_sf::linear::quad(double x, double y, double x1, double x2, double y1, double y2, 
	unsigned int idx, valtype_e vt)
{
	double square;
	switch (10*idx + vt)
	{
		case 0:
			square = (x2 - x1)*(y2 - y1);
			return (x - x2)*(y - y2)/square;
		case 1:
			square = (x2 - x1)*(y2 - y1);
			return (y - y2)/square;
		case 2:
			square = (x2 - x1)*(y2 - y1);
			return (x - x2)/square;
		case 10:
			square = (x2 - x1)*(y2 - y1);
			return (x - x1)*(y2 - y)/square;
		case 11:
			square = (x2 - x1)*(y2 - y1);
			return (y2 - y)/square;
		case 12:
			square = (x2 - x1)*(y2 - y1);
			return -1.0*(x - x1)/square;
		case 20:
			square = (x2 - x1)*(y2 - y1);
			return (x - x1)*(y - y1)/square;
		case 21:
			square = (x2 - x1)*(y2 - y1);
			return (y - y1)/square;
		case 22:
			square = (x2 - x1)*(y2 - y1);
			return (x - x1)/square;
		case 30:
			square = (x2 - x1)*(y2 - y1);
			return (x2 - x)*(y - y1)/square;
		case 31:
			square = (x2 - x1)*(y2 - y1);
			return -1.0*(y - y1)/square;
		case 32:
			square = (x2 - x1)*(y2 - y1);
			return (x2 - x)/square;
		default:
			std::invalid_argument("nfem_sf::linear::quad(...) got an invalid calculating type.");
			return 0.0;
	}
}

bool nfem_sf::linear::triangle_element_is_valid(double x, double y, double x1, 
	double x2, double x3, double y1, double y2, double y3)
{
	// Indice of the vertexes must be ordered anti-clockwise
	// The condition is set to a little bigger than zero to check if the vertexes are on the same line
	if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= NFEM_SF_ZERO)
	{
		throw std::invalid_argument("nfem_sf::linear::triangle_element_is_valid(...) got invalid definition range.");
		return false;
	}

	// The calculating point must be within the triangle
	if ((x1-x)*(y2-y)-(y1-y)*(x2-x) < -1*NFEM_SF_ZERO || 
		(x2-x)*(y3-y)-(y2-y)*(x3-x) < -1*NFEM_SF_ZERO || 
		(x3-x)*(y1-y)-(y3-y)*(x1-x) < -1*NFEM_SF_ZERO)
	{
		throw std::invalid_argument("nfem_sf::linear::triangle_element_is_valid(...) got an invalid coordinate.");
		return false;
	}

	return true;
}

double nfem_sf::linear::triangle(double x, double y, double x1, double x2, double x3, 
	double y1, double y2, double y3, unsigned int idx, valtype_e vt)
{
	double a, b, c, d;
	double cx[3], cy[3];
	cx[0] = x1; cx[1] = x2; cx[2] = x3;
	cy[0] = y1; cy[1] = y2; cy[2] = y3;

	d = (cx[(idx+1)%3] - cx[(idx+2)%3])*(cy[idx] - cy[(idx+1)%3]) - (cx[idx] - cx[(idx+1)%3])*(cy[(idx+1)%3] - cy[(idx+2)%3]);

	if (vt == 0)
	{
		a = (cy[(idx+2)%3] - cy[(idx+1)%3])/d;
		b = (cx[(idx+1)%3] - cx[(idx+2)%3])/d;
		c = 1.0 - cx[idx]*a - cy[idx]*b;
		return a*x+b*y+c;
	}
	else if (vt == 1)
	{
		return (cy[(idx+2)%3] - cy[(idx+1)%3])/d;
	}
	else if (vt == 2)
	{
		return (cx[(idx+1)%3] - cx[(idx+2)%3])/d;
	}
	else
	{
		throw std::invalid_argument("nfem_sf::linear::triangle(...) got an invalid calculating type.");
		return 0.0;
	}
}

bool nfem_sf::linear::cube_element_is_valid(double x, double y, double z, 
	double x1, double x2, double y1, double y2, double z1, double z2)
{
	if (x1 >= x2 || y1 >= y2 || z1 >= z2)
	{
		throw std::invalid_argument("nfem_sf::linear::cube_element_is_valid(...) got invalid definition range.");
		return false;
	}

	if (x < x1 || x > x2 || y < y1 || y > y2 || z < z1 || z > z2)
	{
		throw std::invalid_argument("nfem_sf::linear::cube_element_is_valid(...) got an invalid coordinate.");
		return false;
	}

	return true;
}

double nfem_sf::linear::cube(double x, double y, double z, double x1, double x2, 
	double y1, double y2, double z1, double z2, unsigned int idx, valtype_e vt)
{
	double vol;
	switch (10*idx + vt)
	{
		case 0:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y2 - y)*(z2 - z)/vol;
		case 1:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(y2 - y)*(z2 - z)/vol;
		case 2:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x2 - x)*(z2 - z)/vol;
		case 3:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x2 - x)*(y2 - y)/vol;
		case 10:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y2 - y)*(z2 - z)/vol;
		case 11:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (y2 - y)*(z2 - z)/vol;
		case 12:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x - x1)*(z2 - z)/vol;
		case 13:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x - x1)*(y2 - y)/vol;
		case 20:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y - y1)*(z2 - z)/vol;
		case 21:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (y - y1)*(z2 - z)/vol;
		case 22:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(z2 - z)/vol;
		case 23:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x - x1)*(y - y1)/vol;
		case 30:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y - y1)*(z2 - z)/vol;
		case 31:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(y - y1)*(z2 - z)/vol;
		case 32:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (y - y1)*(z2 - z)/vol;
		case 33:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x2 - x)*(y - y1)/vol;
		case 40:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y2 - y)*(z - z1)/vol;
		case 41:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(y2 - y)*(z - z1)/vol;
		case 42:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x2 - x)*(z - z1)/vol;
		case 43:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y2 - y)/vol;
		case 50:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y2 - y)*(z - z1)/vol;
		case 51:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (y2 - y)*(z - z1)/vol;
		case 52:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(x - x1)*(z - z1)/vol;
		case 53:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y2 - y)/vol;
		case 60:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y - y1)*(z - z1)/vol;
		case 61:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (y - y1)*(z - z1)/vol;
		case 62:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(z - z1)/vol;
		case 63:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x - x1)*(y - y1)/vol;
		case 70:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y - y1)*(z - z1)/vol;
		case 71:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return -1.0*(y - y1)*(z - z1)/vol;
		case 72:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(z - z1)/vol;
		case 73:
			vol = (x2 - x1)*(y2 - y1)*(z2 - z1);
			return (x2 - x)*(y - y1)/vol;
		default:
			throw std::invalid_argument("nfem_sf::linear::cube(...) got an invalid calculating type.");
			return 0;
	}
}

//bool nfem_sf::linear::tetrahedron_element_is_valid(double x, double y, double z, double x1, double x2, double x3, double x4, 
//	double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4)
//{
//	return true;
//}

double nfem_sf::linear::tetrahedron(double x, double y, double z, double x1, double x2, double x3, double x4, 
	double y1, double y2, double y3, double y4, double z1, double z2, double z3, double z4,  
	unsigned int idx, valtype_e vt)
{
	double x12 = x2-x1, y12 = y2-y1, z12 = z2-z1;
	double x13 = x3-x1, y13 = y3-y1, z13 = z3-z1;
	double x14 = x4-x1, y14 = y4-y1, z14 = z4-z1;
	double crx = y12*z13-z12*y13;
	double cry = z12*x13-x12*z13;
	double crz = x12*y13-y12*x13;
	double scl = (x14*crx+y14*cry+z14*crz);
	// No need to get the absolute value here
	//if (scl < 0) scl = -1.0*scl;

	double a, b, c, d;
	switch (10*idx + vt)
	{
		case 0:
			a = -1.0*jacoby3(1.0, y2, z2, 1.0, y3, z3, 1.0, y4, z4)/scl;
			b = jacoby3(1.0, x2, z2, 1.0, x3, z3, 1.0, x4, z4)/scl;
			c = -1.0*jacoby3(1.0, x2, y2, 1.0, x3, y3, 1.0, x4, y4)/scl;
			d = jacoby3(x2, y2, z2, x3, y3, z3, x4, y4, z4)/scl;
			return a*x+b*y+c*z+d;
		case 1:
			return -1.0*jacoby3(1.0, y2, z2, 1.0, y3, z3, 1.0, y4, z4)/scl;
		case 2:
			return jacoby3(1.0, x2, z2, 1.0, x3, z3, 1.0, x4, z4)/scl;
		case 3:
			return -1.0*jacoby3(1.0, x2, y2, 1.0, x3, y3, 1.0, x4, y4)/scl;
		case 10:
			a = jacoby3(1.0, y1, z1, 1.0, y3, z3, 1.0, y4, z4)/scl;
			b = -1.0*jacoby3(1.0, x1, z1, 1.0, x3, z3, 1.0, x4, z4)/scl;
			c = jacoby3(1.0, x1, y1, 1.0, x3, y3, 1.0, x4, y4)/scl;
			d = -1.0*jacoby3(x1, y1, z1, x3, y3, z3, x4, y4, z4)/scl;
			return a*x+b*y+c*z+d;
		case 11:
			return jacoby3(1.0, y1, z1, 1.0, y3, z3, 1.0, y4, z4)/scl;
		case 12:
			return -1.0*jacoby3(1.0, x1, z1, 1.0, x3, z3, 1.0, x4, z4)/scl;
		case 13:
			return jacoby3(1.0, x1, y1, 1.0, x3, y3, 1.0, x4, y4)/scl;
		case 20:
			a = -1.0*jacoby3(1.0, y1, z1, 1.0, y2, z2, 1.0, y4, z4)/scl;
			b = jacoby3(1.0, x1, z1, 1.0, x2, z2, 1.0, x4, z4)/scl;
			c = -1.0*jacoby3(1.0, x1, y1, 1.0, x2, y2, 1.0, x4, y4)/scl;
			d = jacoby3(x1, y1, z1, x2, y2, z2, x4, y4, z4)/scl;
			return a*x+b*y+c*z+d;
		case 21:
			return -1.0*jacoby3(1.0, y1, z1, 1.0, y2, z2, 1.0, y4, z4)/scl;
		case 22:
			return jacoby3(1.0, x1, z1, 1.0, x2, z2, 1.0, x4, z4)/scl;
		case 23:
			return -1.0*jacoby3(1.0, x1, y1, 1.0, x2, y2, 1.0, x4, y4)/scl;
		case 30:
			a = jacoby3(1.0, y1, z1, 1.0, y2, z2, 1.0, y3, z3)/scl;
			b = -1.0*jacoby3(1.0, x1, z1, 1.0, x2, z2, 1.0, x3, z3)/scl;
			c = jacoby3(1.0, x1, y1, 1.0, x2, y2, 1.0, x3, y3)/scl;
			d = -1.0*jacoby3(x1, y1, z1, x2, y2, z2, x3, y3, z3)/scl;
			return a*x+b*y+c*z+d;
		case 31:
			return jacoby3(1.0, y1, z1, 1.0, y2, z2, 1.0, y3, z3)/scl;
		case 32:
			return -1.0*jacoby3(1.0, x1, z1, 1.0, x2, z2, 1.0, x3, z3)/scl;
		case 33:
			return jacoby3(1.0, x1, y1, 1.0, x2, y2, 1.0, x3, y3)/scl;
		default:
			throw std::invalid_argument("nfem_sf::linear::tetrahedron(...) got an invalid calculating type.");
			return 0;
	}
}